#include <UTILS/SBuffer.h>
#include <PECAgent/PEContext.h>
#include <DPS/PH/DpsFunctions.h>
#include "UDOP_OpLocal.h"
#include "letkf_mod.h"
#include "analysis.h"
#include "utilities.h"

SPC_NAMESPACE_USE;
UTILS_NAMESPACE_USE;
using namespace DPS;
using namespace std;
using namespace estd;

/*
 * Function used to initialize the operator
 */
void UDOP_OpLocal::initializeOperator() 
{
  // Get the paramter's value
    char str10[150];
	const DpsOPContext &context = this->getContext();
	const std::string &args = context.getCmdArguments();
	strcpy(str10, args.c_str());
	sscanf(str10,"%*s %*c%d%*c %*s %*c%d%*c %*s %*c%d%*c",&num_ensemble, &radius,&astep);
	
  pxen_bar = (double **)malloc(num_ensemble*sizeof(double *) );
  pxen_mass = (double **)malloc(num_ensemble*sizeof(double *) );
  pyen_con = (double  **)malloc(num_ensemble*sizeof(double *) );
  pxen_tlat = (float **)malloc(num_ensemble*sizeof(float *) );
  pxen_tlon = (float **)malloc(num_ensemble*sizeof(float *) );
  pxen_zlvl = (float **)malloc(num_ensemble*sizeof(float *) );
  pxen_sigh = (float **)malloc(num_ensemble*sizeof(float *) );
  pxen_sigw = (float **)malloc(num_ensemble*sizeof(float *) );
  pxen_sigv = (float **)malloc(num_ensemble*sizeof(float *) );
  pxen_page = (int16_t **)malloc(num_ensemble*sizeof(int16_t *) );
  pxen_hdwp = (int16_t **)malloc(num_ensemble*sizeof(int16_t *) );
  pxen_ptyp = (int16_t **)malloc(num_ensemble*sizeof(int16_t *) );
  pxen_pgrd = (int16_t **)malloc(num_ensemble*sizeof(int16_t *) );
  // Create a user defined metric
  // opm.createMetric("intuples","Number of ingested tuples");
  int i = 0;
  for (i=0; i<num_ensemble; i++)
  {
  pxen_bar[i] = NULL;
  pxen_mass[i] = NULL;
  pyen_con[i] = NULL;
  pxen_tlat[i] = NULL;
  pxen_tlon[i] = NULL;
  pxen_zlvl[i] = NULL;
  pxen_sigh[i] = NULL;
  pxen_sigw[i] = NULL;
  pxen_sigv[i] = NULL;
  pxen_page[i] = NULL;
  pxen_hdwp[i] = NULL;
  pxen_ptyp[i] = NULL;
  pxen_pgrd[i] = NULL;
  }
}

/*
 * Function used to finalize the operator
 */
void UDOP_OpLocal::finalizeOperator() 
{
	int i = 0;
	for (i=0; i<num_ensemble; i++)
	  {
		  if (pxen_mass[i] != NULL)
		  	{
			free(pxen_bar[i]);	
			free(pxen_mass[i]);
			free(pxen_tlat[i]);
			free(pxen_tlon[i]);
			free(pxen_zlvl[i]);
			free(pxen_sigh[i]);
			free(pxen_sigw[i]);
			free(pxen_sigv[i]);
			free(pxen_page[i]);
			free(pxen_hdwp[i]);
			free(pxen_ptyp[i]);
			free(pxen_pgrd[i]);
			free(pyen_con);
			}
	  }
	
  free(pxen_bar);
  free(pxen_mass); free(pxen_tlat);
  free(pxen_tlon); free(pxen_zlvl);
  free(pxen_sigh); free(pxen_sigw);
  free(pxen_sigv); free(pxen_page);
  free(pxen_hdwp); free(pxen_ptyp);
  free(pxen_pgrd); free(pyen_con);
}

/*
 * Function used to process command line arguments
 */
void UDOP_OpLocal::processCmdArgs(const string & /*args*/) 
{
  SPCDBG(L_INFO, "Process command line arguments...", "dpsop");
  // ArgContainer arg(args); // include "ArgContainer.h"
  // arg.argc and arg.argv are now available for use
  // do other initialization here
}


/*
 * Process function for input port 0
 */
void UDOP_OpLocal::processInput0(const IPort0_t & tuple)
{
  SPCDBG(L_DEBUG, "Processing tuple from input port 0 " << tuple, "dpsop");
  int i = 0;
  int j = 0;
  // Obtain values from tuples from streams GlobalX and GlobalY
  dlist<double> xb_bar = (const dlist<double>) (tuple.get_xb_bar() );
  dlist<float> xb_tlat = (const dlist<float>) (tuple.get_xb_tlat() );
  dlist<float> xb_tlon = (const dlist<float>) (tuple.get_xb_tlon() );
  dlist<float> xb_zlvl = (const dlist<float>) (tuple.get_xb_zlvl() );
  dlist<float> xb_sigh = (const dlist<float>) (tuple.get_xb_sigh() );
  dlist<float> xb_sigw = (const dlist<float>) (tuple.get_xb_sigw() );
  dlist<float> xb_sigv = (const dlist<float>) (tuple.get_xb_sigv() );
  dlist<int16_t> xb_page = (const dlist<int16_t>) (tuple.get_xb_page() );
  dlist<int16_t> xb_hdwp = (const dlist<int16_t>) (tuple.get_xb_hdwp() );
  dlist<int16_t> xb_ptyp = (const dlist<int16_t>) (tuple.get_xb_ptyp() );
  dlist<int16_t> xb_pgrd = (const dlist<int16_t>) (tuple.get_xb_pgrd() );
  
  dmatrix<double> en_bar = (const dmatrix<double>) (tuple.get_en_bar());
  dmatrix<double> en_mass = (const dmatrix<double>) (tuple.get_en_mass());
  dmatrix<float> en_tlat = (const dmatrix<float>) (tuple.get_en_tlat() );
  dmatrix<float> en_tlon = (const dmatrix<float>) (tuple.get_en_tlon() );
  dmatrix<float> en_zlvl = (const dmatrix<float>) (tuple.get_en_zlvl() );
  dmatrix<float> en_sigh = (const dmatrix<float>) (tuple.get_en_sigh() );
  dmatrix<float> en_sigw = (const dmatrix<float>) (tuple.get_en_sigw() );
  dmatrix<float> en_sigv = (const dmatrix<float>) (tuple.get_en_sigv() );
  dmatrix<int16_t> en_page = (const dmatrix<int16_t>) (tuple.get_en_page() );
  dmatrix<int16_t> en_hdwp = (const dmatrix<int16_t>) (tuple.get_en_hdwp() );
  dmatrix<int16_t> en_ptyp = (const dmatrix<int16_t>) (tuple.get_en_ptyp() );
  dmatrix<int16_t> en_pgrd = (const dmatrix<int16_t>) (tuple.get_en_pgrd() );

  dlist<double> yb_bar = (const dlist<double>) (tuple.get_yb_bar() );
  dlist<float> yb_clat = (const dlist<float>) (tuple.get_yb_clat() );
  dlist<float> yb_clon = (const dlist<float>) (tuple.get_yb_clon() );
  dmatrix<double> en_con = (const dmatrix<double>) (tuple.get_en_con() );

  dlist<double> y_obs = (const dlist<double>) (tuple.get_y_obs() );
  dlist<float> obs_olat = (const dlist<float>) (tuple.get_obs_lat() );
  dlist<float> obs_olon = (const dlist<float>) (tuple.get_obs_lon() );
  
  //Declare array pointers for function calls. No extra data copy.
  double *px_bar = &xb_bar.front();
  float *px_tlat = &xb_tlat.front();
  float *px_tlon = &xb_tlon.front();
  float *px_zlvl = &xb_zlvl.front();
  float *px_sigh = &xb_sigh.front();
  float *px_sigw = &xb_sigw.front();
  float *px_sigv = &xb_sigv.front();
  int16_t *px_page = &xb_page.front();
  int16_t *px_hdwp = &xb_hdwp.front();
  int16_t *px_ptyp = &xb_ptyp.front();
  int16_t *px_pgrd = &xb_pgrd.front();

  int numpar = xb_bar.size();


  double *py_bar = &yb_bar.front();
  float *py_clat = &yb_clat.front();
  float *py_clon = &yb_clon.front();

  double *obs = &y_obs.front();
  float *obs_lat = &obs_olat.front();
  float *obs_lon = &obs_olon.front();

  //int numpar = xb_bar.size();
  int numline = yb_bar.size();
  int numobs = y_obs.size();
  int16_t tstep = tuple.get_ts();
  bool flag = tuple.get_flag();

    for(i=0; i<num_ensemble; i++)
	{
		pxen_bar[i] = (double *)realloc(pxen_bar[i], numpar*sizeof(double) );
		pxen_mass[i] = (double *)realloc(pxen_mass[i], numpar*sizeof(double) );
		pxen_tlat[i] = (float *)realloc(pxen_tlat[i], numpar*sizeof(float) );
		pxen_tlon[i] = (float *)realloc(pxen_tlon[i], numpar*sizeof(float) );
		pxen_zlvl[i] = (float *)realloc(pxen_zlvl[i], numpar*sizeof(float) );
		pxen_sigh[i] = (float *)realloc(pxen_sigh[i], numpar*sizeof(float) );
		pxen_sigw[i] = (float *)realloc(pxen_sigw[i], numpar*sizeof(float) );
		pxen_sigv[i] = (float *)realloc(pxen_sigv[i], numpar*sizeof(float) );
		pxen_page[i] = (int16_t *)realloc(pxen_page[i], numpar*sizeof(int16_t) );
		pxen_hdwp[i] = (int16_t *)realloc(pxen_hdwp[i], numpar*sizeof(int16_t) );
		pxen_ptyp[i] = (int16_t *)realloc(pxen_ptyp[i], numpar*sizeof(int16_t) );
		pxen_pgrd[i] = (int16_t *)realloc(pxen_pgrd[i], numpar*sizeof(int16_t) );
		pyen_con[i]  = (double *)realloc(pyen_con[i], numline*sizeof(double) );
	}

  if (flag)// && (tstep%astep==0) )
  {
	for (i=0; i<num_ensemble; i++)
		{
			for (j=0; j<numpar; j++)
			{
				pxen_bar[i][j]  = en_bar(i,j);
				pxen_mass[i][j] = en_mass(i,j);
				pxen_tlat[i][j] = en_tlat(i,j);
				pxen_tlon[i][j] = en_tlon(i,j);
				pxen_zlvl[i][j] = en_zlvl(i,j);
				pxen_sigh[i][j] = en_sigh(i,j);
				pxen_sigw[i][j] = en_sigw(i,j);
				pxen_sigv[i][j] = en_sigv(i,j);
				pxen_page[i][j] = en_page(i,j);
				pxen_hdwp[i][j] = en_hdwp(i,j);
				pxen_ptyp[i][j] = en_ptyp(i,j);
				pxen_pgrd[i][j] = en_pgrd(i,j);
			}
			for (j=0; j<numline; j++)
			{
				pyen_con[i][j] = en_con(i,j);
			}
		}


  //Declare variables for submission
  //Call the letkf local analysis function
    letkf_local(num_ensemble, radius, numpar, numline, numobs, px_bar, px_tlat, px_tlon, px_zlvl, py_bar, py_clat, py_clon, obs, obs_lat, obs_lon, pxen_bar, pxen_mass, pyen_con);

    write_parinit(num_ensemble, numpar, px_bar, px_tlat, px_tlon, px_zlvl, px_sigh, px_sigw, px_sigv, px_page, px_hdwp, px_ptyp, px_pgrd, pxen_bar, pxen_tlat, pxen_tlon, pxen_zlvl, pxen_sigh, pxen_sigw, pxen_sigv, pxen_page, pxen_hdwp, pxen_ptyp, pxen_pgrd);
  }

  OPort0_t otuple; 
  otuple.set_ts(tstep);
  otuple.set_flag(flag);
  submit0(otuple);

  // updatePerformanceCounters(opm);

}

/*
 * Punctuate function
 */
/*void UDOP_OpLocal::processPunctuation(const Punctuation::Value & value, unsigned int input)
{
  SPCDBG(L_DEBUG, "Processing punctuation from input port " << input, "dpsop");
  // punctuate0(value);
}
*/

/*
 * Function used to notify the operator about the readiness of all ports (this
 * is useful when submitting tuples from within a context other than the process
 * functions of the operator).
 */ 
/*
void UDOP_OpLocal::allPortsReady() 
{
} 
*/

/*
 * Terminate function, called asynchronously
 */
/*
void UDOP_OpLocal::prepareToTerminate() 
{
}
*/

/*
 * Function used to update performance counters
 */
/*
static void updatePerformanceCounters(DpsOPMetrics & opm) 
{ 
  static long long ntuples = 0;
  // This is not the ideal way of using performance counters
  // Ideally one should store the result from getMetricByName, 
  // in a class member variable. For instance:
  //   DpsMetric & m = opm.getMetricByName("intuples");
  // m could be made a member variable of UDOP_OpLocal
  // Then, the following is a more efficient way to use the counter:
  //   m.setValue(ntuples++);
  opm.getMetricByName("intuples").setValue(ntuples++);
}
*/

/*
 * Checkpointing function for saving state, called asynchronously 
 */
/*
void UDOP_OpLocal::getCheckpoint(SBuffer & checkpoint) 
{
  SPCDBG(L_DEBUG, "Getting operator state...", "dpsop");
}
*/

/*
 * Checkpointing function for restoring state, called asynchronously
 */
/*
void UDOP_OpLocal::restoreCheckpoint(SBuffer & checkpoint) 
{
  SPCDBG(L_DEBUG, "Restoring operator state...", "dpsop");
}
*/


